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^ ; Abstract 

. ■ Non-equilibrium path integral methods for computing quantum free energy differences are ap- 

Cd ' plied to a quantum particle trapped in a harmonic well of uniformly changing strength with the 

^ . purpose of establishing the convergence properties of the work distribution and free energy as the 

i-rt ! number of degrees of freedom M in the regularized path integrals goes to infinity. The work distri- 

^ \ bution is found to converge when M tends to infinity regardless of the switching speed, leading to 

finite results for the free energy difference when the Jarzynski non-equilibrium work relation or the 
Crooks fluctuation relation are used. The nature of the convergence depends on the regularization 
^S| ■ method. For the Fourier method, the convergence of the free energy difference and work distribu- 

^ \ tion go as 1/M, while both quantities converge as 1/M^ when the bead regularization procedure 

f~^ ' is used. The implications of these results to more general systems are discussed. 
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I. INTRODUCTION 

In the preceding paper yj, a non-equilibrium path integral method for computing free en- 
ergy differences based on combining the Jarzynski equality [2|,|3| and the Crooks fluctuation 
relation 4 , |5l] with the path integral formulation of the quantum mechanical partition sum 
p, lZI, la, |9|] was presented. The path integral representation of the canonical partition func- 
tion is based on mapping a quantum system at finite temperature onto a classical system 
with additional degrees of freedom. A non-equilibrium process can be carried out on this 
isomorphic classical system along a well-defined trajectory in fictitious time. The Jarzynski 
and Crooks relations are valid for such a process, but only under the assumption that the 
work distribution converges as the parameter M of the regularization procedure applied to 
the infinite dimensional path integral goes to infinity. 

A regularization procedure is needed because particles are represented as objects with 
an infinite number of degrees of freedom in the path integral formulation. As a result, 
non-equilibrium dynamical processes in this representation can lead to divergences in non- 
physical quantities, such as the average total Hamiltonian of the particle or the work per- 
formed on the system in the fictitious process, which is a central quantity in the Jarzynski 
and Crooks relations. A regularization procedure restricts the number of degrees of freedom 
to a finite number M and results in finite but M-dependent estimators for quantities of 
interest. In the regularized path integral representation, one finds that the expression for 
the work takes the form of the difference between two quantities which diverge as M — > cxo. 
Having an estimator for a physical quantities that take the form of the difference between 
two diverging quantities is not unusual in the context of path integrals [sl], but this does make 
it important to establish the convergence properties of all relevant estimators. Furthermore, 
even when it can be demonstrated that the regularization procedure leads to convergent 
results, the viability of the non-equilibrium path integral method as a means of computing 
quantum free energy differences is strongly dependent on the rate of convergence of the 
regularized path integral to the exact quantum result. Neither the convergence nor the 
rate of convergence was addressed in detail in Ref. [l!, although strong numerical evidence of 
convergence was presented for a quantum particle in a quartic potential. 

In this paper, the rate of convergence of different regularization procedures is examined 
in detail for the special case of a quantum harmonic oscillator. Harmonic systems have 
the advantage of being often amenable to analytical treatment, allowing for closed form 
and exact solutions. Here, the convergence of the regularization procedure is studied in 
three regimes of the non-equilibrium process, i.e., the quasi-static, the finite time and in- 
stantaneous switching. The work distribution is computed for each of these regimes and its 
convergence, as well as the convergence of free energy difference are analyzed as the number 
of degrees of freedom goes to infinity. It will be shown that the free energy difference and the 
work distribution converge for both the Fourier and the bead regularization procedures[l|], 
with the latter converging more quickly. 

The paper is organized as follows: Section [TTl presents a brief overview of the method as 
it applies to the harmonic oscillator. Section IIIII contains the analysis of the convergence 
of the free energy under different regularization schemes. The work distributions will be 
determined using a generating function technique explained in Sec. [IVl In Sections |Vl |Vl] 
and IVIII the quasi-static, finite time and instantaneous switching processes, respectively, 
are studied. A comparison between the non-equilibrium work distribution generated by 
the fictitious dynamics and that generated with real time quantum dynamics is made in 



Sec. IVIIII The conclusions are given in Sec. IIXI 



II. METHOD AND MODEL SYSTEM 



We consider a one- dimensional quantum system with a Hamiltonian operator H{X) = 
T + V, with T = fP' /{2m) and V = V{x, A) = \m\iP'. Here the potential energy V depends 
on a control parameter A, which is equal to the square of the frequency uj. The canonical 
partition function of this system at an inverse temperature (3 is defined by 



Z{X) 



-I3F(\) ^ rj,^ g-/3H(A) 



and can be written asp, \^ 
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where the integral is over closed paths x{s) [i.e., x{[3Ti) = x{0)] and the Euclidean action S 
is a functional of x given by 



S[x, A] 



/3h 



ds 
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Here and below the s dependence of x in integrals over s will always be implied. For this 
one-dimensional harmonic system, the quantum free energy is known to be exactly 



F(A) =(3-Hog[2smh{pnuj/2)]. 



(4) 



The non-equilibrium path integral approach of computing free energy differences uses a 
non-equilibrium process defined by a fictitious dynamics in which A is changed from A^ = i^A 
to Xb = oj\ o^^^ ^ time r, while starting at canonical equilibrium corresponding to A = A^. 
This fictitious dynamics is derived by introducing a new field p{s) which is also periodic in 
imaginary time, satisfying p{s) = p{s + (3h), leading to an expression that has the form of a 
classical partition function: 



Z(A) = C fvxVp e 



■fiH\x,-p,\\ 



where the fictitious Hamiltonian is given by 
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du 
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Here u = s/{(3K) is a scaled imaginary time variable and the string tension is 

m 



2h^- 



j3^h 

We are interested in a Hamiltonian process, with equations of motion 

dx 5H[x,p,X\ p 



dt 5p{u) m 

dp 5H[x,p,\] d'^x 
dt 5x{u) du^ 
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(7a) 
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in which A is time-dependent, and satisfies the boundary conditions A(0) = A^ and A(r) = 
Xb- Then, defining the fictitious work as 

W = H[x{t),p{t),Xb] - H[x,p,Xa], (8) 

the following identities were shown to hold[l| 

{e-n,. = e-'""', (9) 

and 

Pf{W) = e^^e-f'^^Pri-W), (10) 

where (.)aa denotes an average over non-equilibrium trajectories whose initial conditions 
are drawn from a canonical distribution with A = A^, AF = F{Xb) — F{Xa), PfiW) is the 
probability density to do an amount of work W during the process that takes A from A^ to 
Xb (the forward process), and Pr{—W) is the probability density to do work —W during 
a similar process that starts at Xb and ends at Aa (the reverse process). Equations ([9]) 
and fITOl) are the Jarzynski's non-equilibrium work relation and Crooks fiuctuation relation, 
respectively. 

The distribution of the work done as a result of the change in u will be studied for this 
model. This distribution can be used either in the Jarzynski relation (jH]) or in the Crooks 
fiuctuation relation fITU]) to determine the free energy difference, both of which should yield 
AF = Fb~Fa = p-^ log [smh{phu;B/2)/ smh{phuA/2)]. 

The work in Eq. ([H]) is expressed as the difference of the fictitious Hamiltonian at two 
times. The average value of the fictitious Hamiltonian diverges in canonical equilibrium due 
to the infinite number of degrees of freedom in the path integral ([2j), which might pose a 
problem for the very definition of the work distributions. To investigate whether the work 
distribution is well defined, it is necessary to limit the system to a finite number of degrees 
of freedom. In Ref. [l|, two different regularization methods of reducing the path integral 
representations to a finite number of dimensions were introduced. While both schemes are 
similar in their Fourier representations, in one case the degrees of freedom correspond to low- 
frequency modes of the continuous closed string, whereas in the other case they represent 
a discretized lattice version of the string. In the first case the regularization is based on 
statistical arguments motivated by the form of the resulting Hamiltonian, and in the second 
case the regularization is introduced in the lattice representation by means of the Trotter 
formula. 

For finite M, the value of AF found using the work distribution differs from the exact 
quantum result. In fact, for the harmonic oscillator, one can express the free energy explicitly 
as an expansion in inverse powers of M, and thus assess the convergence of AF analytically. 
This will be studied first in Sec. IIIH in which alternative regularizations aimed at improving 
the convergence are also discussed. Note that the convergence of the free energy only requires 
equilibrium considerations. Then, in Sees. |Vl |Vl] and IVIIl the convergence of the non- 
equilibrium work distributions is analyzed for the cases of an infinitely slow switching rate, 
a finite switching rate and an instantaneous switching process, respectively. 

III. THE FREE ENERGY UNDER DIFFERENT REGULARIZATIONS 

Analytical considerations of the harmonic system proceed most easily in the Fourier 
representation. As explained in Ref. [l|, the Fourier transformation takes on slightly different 



forms in the Fourier and the bead regularization methods. In both cases, though, the 
Hamiltonian assumes the form 

H=J2Hk, (11) 

|fc|<fcc 



where kc is a cut-off wave vector, and 



Hk 
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The dispersion relation for the Fourier regularization is given by 



ujk = 27rk\j — 
m 



27Tk 



while that for the bead regularization is 



2M . 7rk 



(12) 
(13) 

(14) 
(15) 



Note that in the Fourier regularization, k^. is a chosen cut-off, whereas from the Fourier 
transform of the bead regularization, we have k^ = \_{M — 1)/2J or, for M odd, M = 2kc + 1. 
Thus, the two regularization methods can be parameterized either by kc or M. 
Given the Hamiltonian in Eq. flT2|) . the partition sum may be written as 



Zm(A) = C Y[ dpkdxk e~ 

\k\<kj 



/3Hfc 



(16) 



where C is independent of A. Each integral is Gaussian and can be explicitly evaluated. 
Since each mode occurs twice in the product (as k and —k) except for A; = 0, one finds for 
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This product must be evaluated separately for the two different dispersion relations of the 
Fourier and the bead regularization schemes in Eqs. ( fT4l) and (ITSll . respectively. 

In the Fourier regularization with the dispersion relation given by Eq. fll4p . Eq. (1171) can 
be evaluated in the limit kc ^ oo hj writing 



-I3AF ^ ^ 
UJB 
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k=l 
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-I fc'=i >- 



' 2'Kk' J 



(18) 



and using the identity [1 



sinh^; 



fc=i 



n 1+ 



f^2^2 



to finally obtain the exact quantum result 



~ smh{l3hujB/2y 



(19) 



It is straightforward to show using Eq. (ITTl) that the limit is approached as k~^ = 0{M~^). 
For the bead regularization, on the other hand, one uses Eq. (1151) and M = 2kc + 1 to 
write Eq. (TT7I) as 



e-^^^= lim 
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One then uses a different identity, namely jl II] 
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(20) 
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to arrive at 



-/3AF 



lim 

Af^oo 



cosh{Marccosh[l + |(^) ]} - 1 
_cosh{Marccosh[l + |(^)^]} - 1 



1/2 



(22) 



which reduces to Eq. ( 11911 as well when the limit M ^ oo is evaluated, with correction terms 
of order M~^ (see e.g. Ref. ill, ). 

There are other regularization methods possible which lead to even faster convergence 
of the free energy at the expense of a more complicated regularized Hamiltonian. Such 
higher-order schemes can be derived systematically by exploiting the analogy between the 
bead representation and the Hamiltonian-splitting method used to obtain integrators in 
molecular dynamics. In fact, the basis of the bead regularization is the splitting form 
in Eq. (30) of the preceding paper [l|, and that same form is also the basis of the Verlet 
scheme to integrate the equations of motions for classical Hamiltonian systems in molecular 
dynamics simulations [l2j. The advantage of using splitting schemes to derive integrators 
for molecular dynamics is that the approximate dynamics is still symplectic, causing such 
integration schemes to be very stable. The analogous property to symplecticity in the path 
integral context is the Hermitian nature of the Boltzmann operator, which is preserved in 
splitting schemes. 

In an attempt to reduce the error due to operator splitting, many alternative splitting 
schemes for molecular dynamics simulations have been derived (see Ref. il3i and references 
therein). Some of these schemes raise the order of the splitting approximation to 0{5t'^) or 
higher. Such splitting schemes can also be used for path integrals and result in a fictitious 
Hamiltonian in which the beads are not all equivalent or which contains explicit correction 
terms of order 0{h ) in the potential, and whose partition function converges to the real 
quantum partition sum as 0{M~'^) or higher. Other splitting schemes derived for molecular 
dynamics simulations are aimed at reducing the error by minimizing the pre-factors in front 



of the leading correction ternis[l4|]. These methods, however, are based on the assumption 
that different correction terms contribute independently and equally to the error. This 
approach has proved useful in molecular dynamics simulations. For instance, it has been 
demonstrated that an optimized second order scheme called H0A2 can often outperform 
the Verlet scheme in molecular dynamics simulations of rigid water molecules 15| . 



We will briefly consider the H0A2 operator-splitting scheme applied as a means to reg- 
ularize imaginary-time path integrals, since it is the simplest splitting-method that may 
improve the rate of convergence of the regularization. It is based on the following operator- 
splitting scheme: 

^-m/M ^ ^-r^|3V/M^-pf/{2M)^-{l-2r1)pV/(2M)^-pf/(2M)^-n|3V/M _^ Q(^M~^). (23) 

For 7] = 1/4, this splitting scheme reduces to the path integral analog of the Verlet scheme of 
molecular dynamics (applied twice), which corresponds to the standard bead regularization 
procedure. However taking a different value of 77, namely 77 ^0.1931833275037836, mini- 
mizes the sum of the squares of the error terms of (9(M~'^)[14]. It is shown in Appendix 
Rl that when applied to the dynamics of a classical harmonic oscillator, this scheme con- 
serves energy better than the Verlet scheme. The application of the H0A2 splitting scheme 
to path integrals is carried out by representing the Boltzmann operator by M factors of 
exp{—PH/M), taking its trace to get the partition function, using the splitting scheme (123]) 
for each factor, inserting completeness relations between each exponential, and perform- 
ing the momentum integration. This procedure leads to a regularized path integral with 
Hamiltonian 



M 



jvi- f ]\/f 1 ^ 

n=l ^ 1"^ ^ 

Here, w\ = Ar], W2 = 2—47] and the Un and Vn are the positions of the odd and the even beads, 
respectively. Odd and even beads are no longer identical in nature when r] 7^ 1/4 because 
the completeness relation inserted either in front of the e"'-^"^''-'^^/'-^*^^^ or the e~^^^^^ in 
Eq. fl2^ gives rise to different potential strengths W2 and Wi. The partition sum Zm = 
Jd'^^ud^'v exp(— /^i^A/) is of Gaussian form, since Hm can be written as Hm = F^VF, 
with F = {ui,Vi,U2,V2, ■ ■ ■)■ The integral can be evaluated in terms of the determinant 
of the matrix V. For the Verlet-like Hamiltonian, a Fourier transform could be used to 
diagonalize the matrix (i.e., to decouple the modes), which facilitates the determination of 
the eigenvalues and thus the determinant of V. Here, a Fourier transform only yields a 
partial diagonalization, i.e., defining 



M 



(where the factor v2 and the shift in the phase in front of f„ are introduced for convenience), 
one gets 

Eu = Y.{ulvt)-'^>.\l] (26) 
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Since each term in the Hamiltonian is simply a 2 x 2 quadratic form, the partition sum can 
be written as a product of two-dimensional Gaussian integrals, each of which is proportional 
to l/i/det Wfc. The finite-M partition sum thus becomes 



- = n 7 —-. —-. TTTIT^- (27) 



Using the identity in Eq. (12T|) gives the result 



1 



coshJMarccosh [l + \{^Y + ^(^) ]} 



(28) 



By expanding Zm in M, one can write this as 



^M 



2 smh{(3huj/2) 



tanh(/3^cu/2)M2 ^ ^ ' 



where the prefactor for the first correction term is z{ri) = 1/6 — ?7 + 21]'^. Thus, different 
values of t] lead to different convergence properties since z{ri) depends explicitly on rj. The 
leading order correction is minimized by choosing z{ri) = 0, i.e., rj = 1/4, which corresponds 
to the standard bead regularization procedure. Other choices for rj lead to larger correction 
terms and hence slower convergence. 

We therefore conclude that using a different splitting scheme can be useful for path 
integrals if the order of the approximation is changed, but that optimized splitting methods 
need not yield any improvement, even if they have been shown to be beneficial in the 
context of molecular dynamics simulations. We will therefore work only with the Verlet-like 
Hamiltonian in the remainder of this paper. 

IV. GENERATING FUNCTION OF THE WORK DISTRIBUTION 

The calculation of the work distribution proceeds most easily by first determining its 
generating function (which coincides with its Fourier transform) 



G{u) = dW e'^'^PiW) (29) 

= (e'"[^(^)-^(°)l)^^, (30) 

where H{t) = if (x(t), p(t), A(t)), which in this case can be written as in Eq. (TTTj) . 
The equations of motion are given by 

dxk Pk .01 N 

(31a) 



dt m 

-mQl{t)xk. (31b) 



dPfc ,..o2 



dt 

From the equations of motion, it is apparent that all Fourier modes evolve independently 
with a time dependent frequency, and that each mode contributes an independent term 

Wk = Hkir) - Hk{0) (32) 
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FIG. 1: The two harmonic potentials -^muP'x^ for which the free energy difference and the work 
distribution in switching from one to the other is studied. In both cases, m = 1, while for Va, 

LO = oja = 1/2 and for Vg, u) = lob = 5/4. 



to the total work W = J2\k\<kc'^^'- Furthermore, the modes are also independent in the 
initial canonical distribution function exp[— /3-ff(x, p, cj^)]. Because the generating function 
of the sum of independent term is the product of the generating functions of the different 
terms, we get for G{u) 

G{u) = II Gkiu) (33) 

\k <kc 

where 

/da;fedpfce™^*-'^^*(°) 



Gkiu) 



(34) 



Below, we will investigate the convergence of the work distribution functions Pf and P^ as 
M is taken to infinity. 



V. QUASI-STATIC PROCESS 



A. Work distribution 



Consider first a process in which u is changed infinitely slowly or quasi-statically from 
oja to ojb- Since each mode k is equivalent to a classical harmonic oscillator with frequency 



r2fc, and ^k changes from ^2^(0) = a/cu^ + uj\ to f2,fc(r) = \/o jI + ujp, one can use that 
Hk{t)/Vtk{t) is an adiabatic invariant for harmonic oscillators [1 6"], to obtain for the work 

Wk = [fifc(r)/fifc(0) - l]/ffc(0). (35) 

Eq. (!M|) then gives for the generating functions of mode k 
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(36) 

(37) 
The inverse Fourier transform 



Note that all 7^ have the same sign as ub — oja-, cf. Eq. ([3. 
of Eq. ([36]) is 

Pk{W) = |7fc|e-^'=^e(7fcVr), 

where B is the Heaviside step function. The inverse Fourier transform of G{u) in Eq. ([33 
is a convolution of these exponential functions, which yields 



(38) 






P(^W) = e{'yoW)J2^kiW)e- 



-YkW 
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where the form of rk{W) depends on whether mode k is degenerate or not. For degenerate 
modes with a^ = a-k (i.e. 1 < k < ["^^^^^J), Ti:{W) is a linear function of W, 



TkiW) 
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while for non-degenerate modes {k 
given by 



0, and /c = ^ if M is even), Tk{W) is a constant 



Tk{W) 



\lk\ 
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(40b) 
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As an example, let the frequency of the harmonic oscillator be switched from uja = 1/2 
to uiB = 5/4, while setting m = h = 1 and /3 = 10 (see Fig. [1]). In Fig. [21 the forward 
and reverse work distributions are plotted, using Eq. (!39|) with M = and M = 15, where, 
for the reverse process, the values of ua and ujb are interchanged. The value of M = 15 
was chosen because the work distribution has then already converged up to about 5%, while 
M = corresponds to the classical process. As explained in the preceding paper [l|, the 
work values at crossing points of the forward and reverse distributions should be equal 
to the free energy difference. In Fig. [21 the classical distributions are seen to cross at 
Wc ~ 0.09, which agrees with the prediction AFdassicai = P~^^og{ujB/ujA) = 0.0916..., 
while the quantum distributions cross at Wc ~ 0.36, which agrees with the prediction 



AK 



quantum 



p-'^\og{smh{hpiUB/'2:)/smh{hPujA/2)) = 0.375... within 4%. Better agree- 



ment for the quantum case is obtained by using higher values of M. 
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FIG. 2: (Color online) Classical and quantum forward {Pj) and reverse {Pr) work distributions for 
a quasi-static switching between the potentials in Fig. [H generated from Eq. ([39]) using the bead 
dispersion relation P^ . with M = (classical case) and M = 15 (quantum case, converged up 
to ~ 5%) The circles labeled C and Q are the crossing points of the classical and the quantum 
distributions, respectively, whose W values should coincide with the classical and quantum free 
energy differences, according to the Crooks fluctuation relation p^ . 



B. Convergence of the w^ork distribution 

Although the distribution functions in Fig. [2] suggest a numerical convergence, one can 
prove analytically that they converge by analyzing the cumulants of the distributions, rather 
than the somewhat cumbersome infinite products and sums in Eqs. (I40ap and (140bp . The 
cumulants Hj of the distribution P{W) follow from the generating function as 



rb-j 



d(m) 



logG(w) 



(41) 



M=0 



which means the generating function can be expressed in terms of cumulants as 



logG(n) = ^-|(i«)^ 



(42) 
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The cumulants are therefore also formally related to the free energy difference by the Jarzyn- 
ski relation Q), i.e., 

1 1 °° 

^F = --log(e-^^) = —logCm = J2^{-(3y-\ (43) 

P P j^^ J- 

Cumulants of independent variables are additive, so that 

k\<kc \k\<kc ^k 



Kj 



where the mode cumulants were determined from Eq. fl5B]) . The convergence of the cumulants 
as fcc — >■ oo can be determined now. From Eq. (1371) . one sees that 

where 

Acj2 = wl - u;^, (46) 

and therefore 

The effect of this asymptotic formula is different in the Fourier and in the bead regularization. 
In the Fourier regularization, uik ~ k, so that Kj ~ 1/k'^K Thus one sees that not only do 
all cumulants in Eq. (jH]) converge as fcc — ^ oo, higher order cumulants converge faster than 
lower orders, i.e. 

Kj{k^) - Kj{oo) = 0{l/kl^-^) = 0{l/M^^-^). (48) 

This means that the shape of the work distribution converges faster than the average. How- 
ever, this property turns out not to be robust. One can show that adding a perturba- 
tive quartic term to the potential causes all cumulants to converge as /c^^ in the Fourier 



regularization [171 



The asymptotic convergence of the cumulants is different in the bead regularization 
scheme, as is evident when the sum in Eq. f H^ is first split up into a sum from k = —k* to 
k* and a sum of modes with k* < \k\ < kc = M/2, and the latter is approximated by an 
integral: 



Kj ~ 



^-<-^)'ff)'ri- 



where Cj is the contribution of the k < k* modes. If k* = 0{hPvAu^), Cj can be shown 
to converge quickly as M — > cxd. Using Eq. flTSl) . one gets 



r\j 



2(j-l)! V 8^^ / J^ sin2^(7rg) 
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FIG. 3: (Color online) Convergence of the quantum forward work distributions for a quasi-static 
switching between the potentials in Fig. [H generated from Eq. (j39p using the bead dispersion 
relation (|15p . Plotted is PjiW) for various values of M. 

where the integration variable was changed to g = k/M. The integral can be performed 



using 



dx' 



sin^-' x' 






(2} - 2)!!(2'i 



cosx 



(2j - l)!!(2z)!! 



sin^*+^ X ' 



leading to 



Kjrj 



Cj 



2(J - 1)! 



(3h^Au^ 



in^ 



2j - 1 A;*2i 



:^ + ^(^" 



(51) 



Thus all cumulants in the bead regularization converge as 1/M^. Note that this behaviour 
is consistent with the numerical results presented in Fig. 3 of the preceding paper [l|]. 

Even though all cumulants converge in the same way, the coefficients in front of the 1/M^ 
terms often turn out to get smaller for higher order cumulants, so that the shape will still 
appear to converge rather quickly, as Fig. [3] illustrates: for large enough M, the shapes of 
the forward work distributions for different values of M are very similar, but shifted along 
the ly axis. 

Since the cumulants of the work distribution converge in the limit M ^ oo in both 
regularization schemes, the work distribution itself is well defined and converges as M~^ or 
M~^ for the Fourier and the bead regularizations, respectively, in spite of the fact that the 
average energy (H) diverges linearly with M. 
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C. Convergence of the Jarzynski relation 

We will end this section with an explicit demonstration that Jarzynski's non-equilibrium 
work relation is obeyed and converges to the correct quantum result in the limit M — »• 
oo. First note from Eq. fl2^ that the left-hand side of Eq. ^ may be reformulated as 
(exp(-/5iy))A^ = G{i(3). Using Eqs. ([MD, dSSD and dSZD, one finds 



Kn 



,2 I , ,2 



-I3w\ = ^ TT ^fc + '^> 



-■■;--n^^- («2) 



This product coincides with the product on the right-hand side of Eq. (fT7|) . showing that 
the Jarzynski equality (e~'^^) = e~^^^ holds and that the resulting free energy difference 
converges to the exact result in the limit M — >■ cx) in the same way, i.e., as 1/M for the 
Fourier regularization and as 1/M^ for the bead regularization, as expected from Eqs. 
and dSI]). 



VI. FINITE SWITCHING SPEED 

A. Work distribution 

Consider now the case in which the switch from oja to ujb is done in a finite (fictitious) 
time r, using the protocol 

cu^(t) =uj\ + Aiu^-, (53) 

T 

where Au"^ was defined in Eq. ( H6l) . so that co'^(t) = cu^. 

As in the quasi-static process, the Fourier modes are independent and can be handled 
separately using a time dependent frequency flk(t) defined through Eqs. (TT3l) and (l53l) . The 
solution of Eqs. (I31al) and (ISlbO for a set of initial conditions X,fc(0) = {xk{0),pk{0)) for this 
switching protocol can be written as a linear mapping 

Xfc(t) = Mfc(t)Xfc(0), (54) 

where Det Mfc(t) = 1 since the dynamics conserves phase space volume. The mapping matrix 



Mfc(t) can be written explicitly in terms of the Airy functions 18|] 0i(t) = Ai(— fi^(t)/6) and 

02(t) = Bi{-Ql{t)/b) with b = (|Acj2|/^)2/3 ^g 

Mfc(t) = Mfc(t)M^i(0) (55) 

where 

''^ ^ \mMt) m(f)2{t)J ^ ' 

which yields 

^ . ^ ( Ai(-i/,) Bi'(-yo) - Bi(-y,) Ai'(-yo) ;;;75[Ai(-y*) Bi(-yo) - Bi{-y,) ki{-y,)]\ 

'^ ' "" \amVh [Bi'(-2/,) Ai'(-yo) - Ai'(-2/,) Bi'(-2/o)] Bi'(-y,) Ai(-yo) - Ai'(-y*) Bi(-i/o) 7 

(57) 
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where cr = ±1 depending on the sign of Acu^, yt = fll{t)/b, and Ai' and Bi' are the derivatives 
of the Airy functions. 

Defining a diagonal matrix 

Imnlit) 



D.(t)=(2 Q 



2m 



the energy contribution of mode k at time t can be written as 

H,it) = X.lit)D,it)X.lit) = Xl{Q)H,{t)Xl{0) (58) 

where the "Hamiltonian" matrix is given by 

H,(t) = M^(t)D,(t)M(t). (59) 

From Eq. flT^ . the work contribution of mode k can then be written as 

W^fc = XnO)[Hfc(r)-Dfc(0)]Xfc(0) 

where we used that Hfc(O) = 0^(0). It is now straightforward to compute the generating 
function Gk{u) of W4, 

r dXfc(O) e^fc (0)[i«Hfc(r)-(/3+m)Dfc(0)]X,(0) 

Gk{u) = - 



/dXfc(0)e-'3x^(o)Dfe(o)Xfc(o) 

-1/2 



Det (l - ^C.) 



(60) 



where I is the identity matrix and 

Q = D,"^/^(0)H,(r)D,"^/^(0)-l. (61) 

To illustrate the dependence of the work distributions on r, we need to perform the Fourier 
inverse on G{u) = Yl^ Gk{u). For \k\ > 0, this is simple, since k and —k are degenerate and 
yield the combined contribution 

^fc(") = T.^w| iuc) = ^ ^' (62) 

Det(l-^Qj i_,^^i_,^^ 

where A^ and A^ are the two eigenvalues of Cfc. The right hand side of Eq. fl62|) has 
the form of two independent exponential modes with rates /3/A^ ' . The Fourier inverse of 
^o<\k\<k,Gk{u) is therefore 



AK(\ 



P{W) = Q{^',W) Y^ ^,k 7^ ^"'^"^' (63) 



Ig = i V -7' 



,W 



where by definition A^ = /5/72(A;-i)+i- This expression does not contain the contribution of 
the mode k = 0, for which 

Goiu) = 



Det (I - |Co) V(l-m/7+)(l-W7- 
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instantaneous forward 
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FIG. 4: (Color online) Forward and reverse quantum work distributions for various switching 
speeds between the harmonic potentials in Fig.[Tl generated from Eq. ([39]) for the quasi-static case, 
from Eq. (j65p for the finite r = 3/2, and from Eq. (j73p for the instantaneous case. Note that the 
intersection points occur at the same value of VF, as they must according to the Crooks fluctuation 
relation. 



.(1) 



.(2) 



where 7+ = I3/\q and 7_ = /?/Aq , and whose Fourier inverse is given by 



Po{W) = v/7T7^e- 



7+ +7- 



W T 

-'0 



7+ + 7- 



W 



(64) 



where Iq is the zeroth order modified Bessel function of the first kind. The derivation of 
this result, which holds when 7_|_ and 7_ are both positive, is given in Appendix [Bl In fact 
it can be shown that for monotonically increasing u, such as in Eq. (1531) . 7+ and 7_ must 
always be positive, though this restriction need not hold for more complicated protocols. 
The complete work distribution function is therefore 



P{W) 



w 



dW Pq{W')P{W -W), 



(65) 



which unfortunately does not lead to a closed form, but can easily be computed numerically 
since Pq and P are known analytically. As an example, the resulting P{W) for the same 
parameters as for the quasi-static case is shown in Fig. IHfor r = 3/2, both for the forward 
and reverse process. 
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B. Convergence of the w^ork distribution 

To determine the convergence properties of the work distribution as M ^ oo, we will 
once again use the cumulants of the work. From the generating function in Eq. (I60!l . we 
have 

g(^.f^4Ti..og(,-|Q). (66) 

By expanding the logarithm and equating like powers of u for j > 1, one finds the cumulants 
of the different modes to be 

'^r = ^^^Trq. (67) 



One thus has for the cumulants of the total work 

|fc|<fc, 



K 



ETrCi 



This is of similar form to the expression (144|) for the work cumulants in the quasi-static 
process, which could be shown to converge because of the asymptotic property of 7^ as 
A; ^ 00 in Eq. (145|) . A similar asymptotic property here has to involve the two eigenvalues 
A^, and A^ of the real symmetric 2x2 matrix C^ for large k. This asymptotic analysis 
starts with the behavior of the Airy functions for large negative argument s|18|] 

Ai'(-y) ~ -f2ilM!±ilj,./4 (esc) 

B.'(-y) ~ 5!liMl±lly'/4_ (68<j) 



with corrections of relative 0{y ^^'^). Using these asymptotic expressions, one can write for 
the mapping matrix needed in Eq. fl6Tl) 

hA f \ ( ~ *"°^ ^^ mU (0) ^^^ ^'^^ r{in\ 

\-r]kmnk{0) sin 9k VkCosOk /' 

with Ok = livr^"^ - 2/0^^) and r]k = ^^fc(r)/fifc(0). Substituting this Mfc(r) into Eq. ([59]) for 
Hfc using r2fc(r) = rilQk{0) in Dfc(r), and substituting the result into expression ( 16T1) for Ck, 
yields 

Cfc~(r]^-l)l. (70) 

Thus it is clear that asymptotically, both eigenvalues of C^ are equal to A^ = ■^i — Vk~^- 
But ril — 1 is precisely equal to P/'jk, so that Eq. (!67l) gives 

(k) U-i)\ f/3y (j-i)! 
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which precisely coincides with the asymptotic expression of the cumulants of mode k in the 
quasi-static case, Eq. (14^ . In retrospect, this is not too surprising: higher k modes have 
higher frequencies Qk, so that the relative switching speed l/(rf2fc) decreases for increasing 
k, resulting in a quasi-static behavior for the work distributions of the high-/c modes. 

Since the behavior of the work cumulants for large k is the same as in the quasi-static 
switching, the cumulants themselves converge for the finite switching speed as well, and in 
exactly the same way, i.e., (9(l/M^-'~^) in the Fourier regularization and (9(1/M^) for the 
bead regularization. However it should be stressed that the correspondence between the 
quasi-static and finite-switching case only holds for the large k modes, and that the lower k 
modes do differ between the two cases and the total work distribution can vary substantially 
depending on the switching speed. 

C. Convergence of the Jarzynski relation 

With these results for the generating function, the Jarzynski equality i^ can be once 
more checked. Setting u = ij3 in Eq. flM?]) . the generating function gives 

■DetDfc(0)V/^ /DetDfc(O)^^/^ 



G^m 



DetHfc(r)y V^et Dfc(r) 



(71) 



fi,(r)' 

where the fact that Det M (r) = 1 has been used. Putting together all the contributions to 
G{iP), one obtains the same expression as for the quasi-static case, Eq. ( I52l) . The subsequent 
calculation therefore applies, showing that the Jarzynski relation, Eq. (I9l), in the form of 
Eq. (IT9l) . also holds for finite switching rates. 

The result in Eq. (TT^ holds not only for any r using the protocol in Eq. (135]) . but in fact 
for any protocol. A different Hamiltonian dynamics would only change the linear mapping 
Mfc(r), which would still obey Det Mfc(r) = 1, since this follows purely from the Hamiltonian 
nature of the dynamics. Thus, Eq. fl7T|) would still hold, from which the Jarzynski relation 
follows. 

VII. INSTANTANEOUS SWITCHING 

While the work distribution for finite r could only be expressed analytically up to a 
convolution, for r = 0, it is possible to derive a fully analytical form. In this instantaneous 
limit, the system does not have time to change its positions or momenta, whence Wk = 
ImAtu^lxfcp. This allows the generating function Gkiu) to be computed, leading to 



G{u) = , ^ ft ^—-^. (72) 

\/ -L - ^u-T-w k=i pni 



i3oj^ «.•=! pn^^ 



The generating function can be inverted, to yield 



' ^ / '"V, ^ -I — r Zf„ / ^ 'V, \ 



fc=l 






18 



where 7^' = f3Ql/Auj'^ and erfi(x) = — ierf(ix) is the complex error function. This work 
distribution P{W) has also been plotted in Fig. HJ for the same parameters as for the 
quasi-static and finite-r cases and both for the forward and reverse process. One sees that 
regardless of the speed of the process, the crossing points Wc of the forward and reverse 
distributions are all the same, and equal to the free energy difference. 

To check the convergence of the Jarzynski equality in the instantaneous switching case, 
one substitutes u = i[3 into Eq. (1721) . For each of the modes, the generating function at 
u = ij3 coincides with the result for the finite switching speed, i.e. the right-hand side of 
Eq. (1711) . As a result, the convergence of the free energy as computed from the Jarzynski 
equality in the instantaneous case is the same as it was for the finite switching case, as would 
be expected since the former is the r — > limit of the latter. 

The convergence of the work distribution might be expected to follow similarly by taking 
the limit r ^ of the cumulants found in the case of finite switching, but it turns out that 
the limits r — > and fc ^ 00 do not commute. The physical reason is that no matter how 
fast the finite switching is, for fixed positive definite r and growing k, there are always k 
modes whose frequencies fik sue faster than the switching rate r~^, and for those modes 
the switching is nearly quasi-static. But if r = 0, then the process is instantaneous for all 
modes, regardless of their k value. 

To see the noncommutativity of the limits mathematically, note that for instantaneous 
process in which the positions and momenta of the system do not have time to change, 
the Mfc(r) matrix is equal to the identity matrix. The matrix Cfe defined in Eq. ( 16T]) then 
assumes the form 

^/ y , (74) 

where fi^^ = ljI + uj\ and Ql^ = lim,-^o ^li'^) = "^fc + ^b- This form is not the same as the 
r — i> limit of Eq. (!70l) . From Eqs. fl67j) and (!74l) . one finds for the mode cumulants in the 
instantaneous switching case 

M _ U - 1)! f^ls 






(j - 1)! f^ 



27i \^ 



^kB 



l-kA 



l]\ (75) 



where Eq. ( 137|) was used. For large k, Vtks/^kA + 1 — > 2, so that k^- oc 7^-' in this 
regime. Apart from a k independent prefactor, the asymptotic behavior of the cumulants 
as a function of k is the same as in the previous cases, so that the convergence of the work 
distribution is also the same. 



VIII. COMPARISON WITH REAL TIME QUANTUM DYNAMICS 

The viability of the non-equilibrium path-integral approach for the computation of quan- 
tum free energy differences is made possible by the use of a fictitious dynamic, thus avoid- 
ing the problems associated with the description of the real time quantum evolution of 
the system, in contrast to other extensions of the Jarzynski and Crooks fluctuation rela- 
tions to the realm of quantum mechanics, based on the operator formulation of quantum 
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FIG. 5: (Color online) Comparison of the work distributions in the non-equilibrium process of 
switching the oscillator strength linearly from loa = 1/2 to ujb = 5/4 in a time r = 3/2, starting 
in canonical equilibrium with inverse temperature /3 = 4, for quantum, classical, and and fictitious 
dynamics. 



dynamics |19l. |20| . |21| . |22| . |23| . |2J, |25|. As a consequence, however, the work computed within 
this scheme bears no relation to the work performed in any real quantum process except 
in the classical limit hf3 ^ 0, where there are no contributions to the work distribution 
from the non-zero k modes in the fictitious dynamics. The work distribution for an isolated 
quantum system in which the frequency is changed in real (rather than fictitious) time has 
recently been worked out by Deffner and Lutz|25], so that a direct comparison between the 
fictitious dynamics and real quantum dynamics is possible. 

The work distribution for the quantum harmonic oscillator consists of a sum of delta 
functions, since the real quantum work distribution can be obtained by summing over all 
possible transitions with the appropriate transition amplitudes worked out by Husimi 261 . 
Thus, for the purpose of comparison, it is better to consider the cumulative distribution func- 
tion (CDF) C{W) = J_^dw P{w)^^. Figure [5] shows the comparison between the classical 
CDF [cf. Eq. flMl) ]. the quantum CDF (using Husimi's transition amplitudes) and the CDF 
resulting from the fictitious path dynamics [cf. Eqs. fl63l) - fl65|) ] for uja = 1/2, ub = 5/4, 
r = 3/2 and (3 = 4. The three curves are obviously very different. The quantum distribu- 
tion is composed of steps, while the other two are continuous. The quantum distribution 
also does not agree with the fictitious dynamics on average. Interestingly, the deviation of 
the real quantum dynamics from the classical case is less then the deviation from the ficti- 
tious dynamics. Another difference between the real quantum dynamics and the dynamics 
in the other cases is the non-zero probability of negative work values predicted by the dis- 



20 



tribution, which, as suggested in Ref. l25|, is due to excited states decaying to lower states. 



This negative tail disappears in the classical limit for the monotonic protocol in Eq. (!53l) . 



The negative tails in the classical limit found in Ref. l25| can only exist for non-monotonic 
protocols. 

Given these observations, it is clear that the fictitious dynamics is very different from 
real quantum dynamics and has little direct physical content. 

IX. CONCLUSIONS 

The non-equilibrium methods for the calculation of free energy differences in quantum 
systems in the context of the path integral representation of the canonical partition function 
presented in the previous paper were applied to the harmonic oscillator to show that the work 
distribution function is well-defined as the regularization parameter M is taken to infinity. 
Instead of using the real quantum dynamics of the system, the path integral representation 
allows a fictitious path to be defined for which the Jarzynski and Crooks relations are valid. 
By evolving the ring polymer in the path integral representation under fictitious dynamics, 
the difficulties associated with the complexity of the full evolution of a quantum system are 
avoided. 

In particular, expressions for the distribution P{W) of the work W done in the non- 
equilibrium fictitious process in which the strength changes linearly during a time r, were 
derived, in the form of a single convolution for finite r, and in fully explicit form for r ^ 
(the instantaneous limit) and r — ;> oo (the quasi-static limit). From P{W), it was shown 
that the Jarzynski relation holds for this case for any dynamic switching process based on 
(isolated) Hamiltonian dynamics. The convergence of the resulting free energy difference 
was obtained for both regularizations, and goes as 0{M~^) for the Fourier regularization 
and as 0{M~'^) for the bead regularization. The nature of the convergence of the cumulants 
of P{W) as M — i> oo was also determined for any r. Whereas the beads regularization leads 
to cumulants which all converge as M~^, the jth cumulant converges as 0{M^~'^^) in the 
Fourier representation, implying that the shape of the distribution converges faster than its 
position along the W axis. However, one can use perturbative arguments with the harmonic 
oscillator as the zeroth order system to show that the work distributions in the Fourier and 
bead regularization converge as 1/M and 1/M^, respectively. Indeed, in the path integral 
simulations using the bead regularization presented in the preceding paper |l|, one sees a 
1/M^ convergence for the free energy difference as well as for the first and second cumulant 
of the work distribution function. Given the analytical proof given in this paper and the 
numerical evidence in the preceding paper, it can be expected that the convergence of the 
method is general. 

Other regularization schemes based on the splitting method were also briefly considered 
and it was found that splitting schemes optimized for molecular dynamics need not be 
optimal for the convergence of path integrals, in contrast to higher order splitting schemes 
which will have better asymptotic convergence properties. 

The difference between the fictitious dynamics and real quantum dynamics was demon- 
strated by direct comparison of the work distribution. Not only is the nature of the dis- 
tribution different (delta peaks for the quantum case, a smooth function for the fictitious 
dynamics), but the two also do not agree on average. Nonetheless, the free energy found 
from the non-equilibrium method with fictitious dynamics is the exact quantum free en- 
ergy difference. 
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APPENDIX A: ENERGY CONSERVATION OF THE CLASSICAL HARMONIC 
OSCILLATOR UNDER THE HOA2 INTEGRATION SCHEME 

The H0A2 scheme is a second order integrator for molecular dynamics simulations, like 
the Verlet scheme, but contains a parameter 77 which can be tuned to make the dynamics 



as accurate as possible |13l. 1 141 . |28|. In this appendix, the level of energy conservation as a 
function of f] will be determined explicitly for the classical harmonic oscillator. 

The state of the classical harmonic oscillator is represented by a position x and momentum 
p. Combining these quantities into a two-dimensional vector F = {Ljx,p/m), the energy of 
the oscillator is given by E = ylTp. 

In applying the H0A2 splitting scheme in Eq. (|23l) to molecular dynamics, the kinetic 
operator T is to be replaced by the free Liouville operator Ct = {p'^/2m, .} and the potential 
operator V by the interaction Liouville operator Cy = {V{f), ■}, where {-, ■} is the Poisson 
bracket operator. The sum of these two Liouville operators is the full Liouvillian C = 
Ct + Cv- Because of the linearity of the equations of motion of the classical harmonics 
oscillator, the exponentials of these two Liouville operators can be written as the matrices 

e^-* = A(t) = (\ :^ (Al) 



^0 1^ 
e^-* = B(t) = (_^^ J) , (A2) 

respectively, which act on the vector F, and where r = ut. For the one-step propagator, one 
thus finds from Eq. (!23|) 



\m) \2m) \ m J \: 



e-'"^Bl^)M^)Br-±^^)M^)B(!l 



^ 2W^, 4M3 M 4M3 \ { K'W 

r_ , »?(l->?)r^ _ >?^(l-2>?)r5 i _ ^^ , r,(l-2>;)r4 I • l^-^J 

M "•" JV/S 4M5 ^ 2M2 "I" 4M4 / 

We will denote this approximate propagator-matrix by P(t/M). The approximate propa- 
gator P(t) over a time t is given by the Mth power of P(t/M), which can be evaluated by 
diagonalization: 

P(^)=U-diag(/ii,/i2)-U-\ (A4) 

with /ii and [I2 the eigenvalues of P(t/M), so that 

P(t) = U-diag(/xf,/xf)-U-^ (A5) 
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Because we want to minimize the error in front of the leading correction term, we only need 
the expansions of the eigenvalues and of U in inverse powers of M. The eigenvalues may be 
expressed as 



/ii = l + 



IT r^ l{l + 27] - A7f)T^ / 1 \ 



1^2 



M 2M2 8M3 

IT r^ i(l + 2?7 - Ar]^)T^ 






M 2M2 



+ 



whence 



/^f = e- 



/if 



^M3 



+ 






1 + 



i(l - 6r/ + 12r/^)r3 , ^^/ 1 \ 



M2 



+ 



VM37 



_ i(l - 6r/ + 127^2)^3 ^ 

M2 VMsy 



while the matrix U can be written as 
U = 



_. _|_ i(l-6q+4??^)T^ ■ _ i(l-6??+4t;^)T ^ 



8M2 



67?+V)t2 \ / 1 \ 



Substituting these expressions into Eq. (jASp . the deviations in the total energy from its 
initial value as a function of the initial conditions xo and po can be written as 



E{t) - E(0) 



lm\Pit)-W 



^ IF |2 

-mi n 
2 ' °' 

6r7 + 4?72)r2 sin^r , ,, , 



^M2 



pl/m — 2iJXoPo cot r) + 0{M 



(A6) 



Thus, the leading violation of energy conservation is of order l/M^ for general rj. In par- 
ticular, the Verlet scheme, for which rj = 1/4, has second order violations in the energy 
conservation. Equation flA6p shows that the leading order violation can be eliminated alto- 
gether by taking the solution of 1 — Gr^ + At]'^ = 0, i.e., r/ = (3 ± \/5)/4. The larger of these 
solutions leads to negative time steps in the H0A2 scheme, which can lead to instabilities, 
so the value of rj to take is 77 = (3 - \/5)/4 = 0.1909830056250526 . . .. This value of 77 
makes the scheme pseudo-fourth order for the harmonic oscillator, and is very close to the 
general optimized value of 77 = 0.1931833275037836 . . .[28], which has been demonstrated to 
have smaller variations in the total energy than the Verlet scheme for general potentials in 
numerical simulations [li 



APPENDIX B: FOURIER INVERSE FOR THE CLASSICAL WORK DISTRI- 
BUTION 



To arrive at the classical result Eq. ( |64l) for the finite-r switching process, one needs to 
compute 



1 r°° 

Po{W) = — du 
2n .1^ 



-iuW 



v/(l-iM/7+)(l-iM/7-) 



(Bl) 
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FIG. 6: Integration contours (dotted lines) and branch cuts of Gq{u) (wiggly lines) for different 
signs of 7+ and 7_ , for the three cases that need to be distinguished in the evaluation of the Fourier 
inverse in the appendix. The crosses are the singular branch points. 



The integrand is a two- valued complex function with two singularities — i7±, which are also 
the branch points. The branch to take in the integral should be that for which the integrand 
is 1 at M = [since Go(0) = / dW^Po(W^) = !]• In addition, no branch cut should be crossed 
as one integrates from — oo to +00, i.e., the branch cut should not cross the real axis. This 
restricts the choice of where to put the branch cut: if 7+ and 7_ are of opposite sign, with 
7+ > say, then in order to avoid a branch cut on the real axis the branch cut should have 
two parts, one extending from —27+ to — ioo and one from i|7_| to zoo, as in the left panel 
of Fig. O If, on the other hand, 7+ and 7_ are both positive, then the branch points — i7± 
both lie on the negative imaginary axis, and it is convenient to put the branch cut between 
the two branch points, as indicated in the middle panel of Fig. [6l while if 7+ and 7_ are 
both negative, the situation is mirrored with respect to the real axis, as shown in the right 
panel of Fig. [61 This third case can also be treated by setting W —^ —W and u — > —u in 
Eq. (IB1|) . and therefore does not need to be treated separately. 

For the first case, i.e., 7+ and 7_ of opposite sign, with 7_ chosen negative, one can shift 
the integration line up or down without hitting the branch points, until the integration line 
passes straight through the middle of the two branch points. In the integral in Eq. fIBip . this 
corresponds to a shift of the integration variable over i(7_|_ + 7-)/2. Performing in addition 
a scaling, such that u = [— i(7+ + 7-) + t|7+ — 7_|]/2, gives 
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This result was 



2 

using formula 8.432.5 in Gradshteyn and Ryzhik with u = and 2: = 1 291 . 
also found by Deffner and Lutz [Eq. (27) in Ref. l25|] when they considered the classical limit 
of a real quantum dynamical process. However, 7+ and 7_ can only have opposite signs if 
the oscillator frequency is changed non-monotonically, which is not the case considered in 
the text. 

If 7+ and 7_ are both positive, which occurs for monotonically increasing oscillator fre- 
quencies [cf. Eq. (l53p ]. one cannot simply shift the integration line to run in between the 
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branch points since one of the branch points would be crossed. The integral is then evalu- 
ated as follows. For W < 0, one can construct a contour Ci composed of the real axis and 
an infinite semi-circle in the upper half of the complex plane, as shown in Fig. [61 Because 
for W < 0, the factor e~'"^ decays exponentially when u has a large positive imaginary 
component, the semi-circle does not contribute to the integral, and the contour integral over 
Ci has the same value as the original integral. Since no singularities lie within this contour, 
the integral is zero, proving that for W < 0, the integral is zero. If, on the other hand, W 
is positive, then the integrand decays for u with a large negative imaginary component, and 
the integral can be replaced by an integration over the closed contour C2 found by adding 
an infinite semi-circle in the lower half of the complex plane, cf. Fig [61 Note that because of 
the location of the branch cut, one remains on the same branch of the function going along 
C2, which is a requirement for the contour to be truly closed. The contour C2 does contain 
singularities, and in particular, it contains the branch cut. The contour can be deformed 
without crossing singularities to the contour C3 in Fig. [^1 which goes around the branch cut. 
One easily shows that the contributions from the parts that go around the branch points 
vanish. Furthermore, the integrand changes sign across the branch cut and the contour is 
traversed in opposite directions on either sides, so that the contributions from the left and 
the right segments of C3 are equal, and the integral becomes 

/-i7- p-iuW 

-i7+ 7rV(l-m/7+)(l-iM/7-) 

Changing integration variables to t where u = [— i(7+ + 7-) + 1^(7+ — 7-)]/2 reduces this 
integral to 

Given the representation of the modified Bessel function of the first kind given by formula 
8.431.1 in Gradshteyn and Ryzhik with i> = 0^], one finds Eq. 
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